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Comparative study between TV-body and Fokker-Planck simulations 
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ABSTRACT 

To understand the effects of the initial rotation on the evolution of the tidally limited clusters 
with mass spectrum, we have performed TV-body simulations of the clusters with different 
initial rotations and compared the results with those of the Fokker-Planck (FP) simulations. 
We confirmed that the cluster evolution is accelerated by not only the initial rotation but also 
the mass spectrum. For the slowly rotating models, the time evolutions of mass, energy and 
angular momentum show good agreements between TV-body and FP simulations. On the other 
hand, for the rapidly rotating models, there are significant differences between these two ap- 
proaches at the early stage of the evolutions because of the development of bar instability 
in TV-body simulations. The shape of the cluster for TV-body simulations becomes tri-axial 
or even prolate, which cannot be produced by the 2-dimensional FP simulations. The total 
angular momentum and the total mass of the cluster decrease rapidly while bar-like structure 
persists. After the rotational energy becomes smaller than the critical value for the bar instabil- 
ity, the shape of the cluster becomes nearly axisymmetric again, and follows the evolutionary 
track predicted by the FP equation. We have confirmed again that the energy equipartiton is 
not completely achieved when M2/ Mi(rri2/ mi) 3 / 2 > 0.16. By examining the angular mo- 
mentum at each mass component, we found that the exchange of angular momentum between 
different mass components occurs, similar to the energy exchange leading to the equipartition. 

Key words: stellar dynamics, globular clusters: general, methods: N-body simulations, meth- 
ods: statistical 



1 INTRODUCTION 

The effects of initial rotation on the dynamical evolution of star 
clusters received substantial attention because the rotation is a nat- 
ural consequence during the formation process. The current popu- 
lation of star clusters may not show significant amount of rotation, 
but it does not mean that the initial conditions inhibit the presence 
of rotation. The direct measurement of the rotation is difficult since 
it requires long integration with large aperture telescope. Rather in- 
direct evidence for the rotation comes from the shape of globular 
clusters. Although most globular clusters show high degree of cir- 
cular symmetry, the el lipticity has been measured for large num ber 
of star clusters (e.g.. lWhite & Shawlll987tlChen & CherfcOlCl) . If 
the ellipticity is due to the rotational flattening, many clusters still 
have some degrees of rotation. Even though the amount of rota- 
tion in current population of globular clusters is rather small, initial 



* E-mail: jshong@astro.snu.ac.kr (JH); eunhyeuk@kari.re.kr (EK); hm- 
lee@snu.ac.ki (HML); spurzem@bao.ac.cn (RS) 



clusters could have been rotating much more rapidly since the ro- 
tation phases out as the clusters evolve dynamically. 

The effects of rotation on t he dynamical ev olution have been 
studied by a number of authors. lGoodmanl(ll983l) has extended the 
Fokker-Planck (FP) equation for rotating systems, but his study was 
limited to slowly rotating systems by imposing the spherical sym- 
metry for the shape of t he clusters. The therm odynamical analyses 
have been pioneered bv lHachisul ( fT97^.ll982h and found that there 
exists an instability similar to gravothermal catastrophe and they 
named this phenomenon as 'gravo-gyro catastrophe' . Theses earlier 
studies provided the basis of the possible acceleration of dynamical 
evolution due to the initial rotation. 

More careful studies for the rotatin g systems have been car- 
ried out bv lKim et af] f2002L [20041. |2008h usi ng 2-Dimensional FP 
code developed bv lEinsel & Spurzeml (1999). These papers inves- 
tigated both isolated and tidally limited clusters, and single and 
multi-component clusters. The general result emerged from these 
studies is that the rotating clusters undergo faster evolution than 
non-rotating ones for single component models. The acceleration 
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is also expected in multi-mass models as well, but the degree of ac- 
celeration could be significantly reduced since the energy exchange 
between different mass components is another accelerating process 
and these two processes compete each other. 

The suitability of the FP approach to the study of dynami- 
cal evolution of star clusters has been a matter of debate because 
the absence of the accurate knowledge on the third integral for ro- 
tating systems prohibits us to include all the possible integrals in 
constructing the FP equation. Comparison with TV-body calculation 
should provide a clue to the validity of the current version of the FP 

i 1 — — in i 1 k 

approach (e.g., Giersz & Heggie 1994a. b; Giersz& Spurzem 1994). 
Such a co mparison fo r rotati ng systems was done by Ernst et al.l 
d2007l) and I Kim et ail d2008l ) for single component models and 
showed that the FP results are generally consistent with the TV-body 
calculations. 

We extend the comparison to the two-component models as an 
interim step to the full multi-mass models. Two-component models 
have the ingredients for the multi-mass models, but have smaller 
model parameters. The important difference between single and 
multi-mass models is the existence of energy exchange among dif- 
ferent mass components. In rotating models, there is also a pos- 
sibility of exchange of the specific angular momentum through the 
dynamical process. It is much easier to investigate such processes in 
TV-body. Another motivation for carrying out TV-body simulations 
and comparing with FP results is to address the validity of the ax- 
isymmetric assumption which is inevitable for the FP approaches. 
If the cluster rotates rather rapidly, the bar-like structure can form 
even with the initial assumption of the axisymmetric shape. The 
evolution of the elongated cluster could be different from the per- 
fectly axisymmetric one. 

This paper is organized as follows: In Section 2, we describe 
the models and their parameters in detail. The effect of the initial 
rotation on the cluster shape will be presented in Section 3. In Sec- 
tion 4, we will compare TV-body and FP results in various angles. 
We will discuss the effects of mass spectra on the dynamical evo- 
lution of star clusters in Section 5. 



2 THE MODELS 

Most of FP results for rotating stellar system with initial mass spec- 
trum are based on the 2D FP solver, m FOPAX dKim et alJ[2004h 
which is the rev ised version of FOPAX dEinsel & Spurzemlll999t 
I Kim et al]| 20021 ) suitable for the rotating stellar system with initial 
mass spectrum. For the complet e description of mFOPAX readers 
are referred sections 2 and 3 of iKim et al.l 112004). 

The NBODY code which we used for this study is one of se- 
ries of direct TV-body programs developed by S. Aarseth since 
1960s. Each version of the NBODY codes has been added some 
epochal schemes such as the Ahmad-Cohen neighbor scheme, 
the K ustaaheimo-Stiefel or chain regularization scheme dAarsethl 
1999). More recently, the NBODY4 and the NBODY6 codes can per- 
form mo re precise calcu lations thanks to the 4th-order Hermite in- 
tegrator dAarsethlll999l) . Specially, the NBODY4 code is designed 
for running on the GRAPE, which is a special-purpose machine for 
only direct TV-body simulations by calculating gravity or coulomb 
interaction with high parallelization. We used the GRAPE6-BLX64 
boards for the TV-body calculations reported in this paper. 

To prepare the initial models for the present TV-body runs we 
have assigned the positions and the velocities of the stars from 
the predefined density, potential and distribution function of ini- 
tial models. These initial models also used in FP runs are ob- 



tained follo wing iLupton & Gunnl Jl987|) and also applied in previ - 
ous studies dEinsel & Spurzem! 1999l;lKim et alj200a.l2004l2008h . 
Globular clusters or open clusters which are typical stellar sys- 
tems considered in the present study are tidally limited by their 
host galaxy. To investigate the time evolution of rotating stellar 
system under the influence of the tidal effect of the host galaxy 
and to compare the dynamical evolution of rotational stellar sys- 
tem using two different numerical methods (FP and TV-body), we 
assume that the stellar systems are orbiting around the center of 
the Galaxy with circular orbit, whic h is already applied in previ- 
ous FP runs dKim et al]|2002l . [20041) . In order to directly compare 
TV-body results with FP methods, we remove the stars whose total 
energy exceeding the tidal energ y due to the host ga l axy instan- 
taneously (i.e., energy cut-o ff, Takahashi et al. |l997t iBaum gardt 
l200lUKidl200l ; lEmst etai]|2007l : lKim et alfeOOgfT According to 
the TV-body computation, it takes at least a crossing time to es- 
cape from the clusters for stars with total energy larger than the 
tidal energy. Therefore, instantaneous removal of stars with en- 
ergy greater than tidal energy is somewhat unrealistic. The prob- 
lem becomes more serious for small- TV systems since the fraction 
of stars to be unbounded at a given time is higher than large- TV 
system s {Takahashi & Portegies Zwarj[T998l ; iFukushige & H eggie 
2000). However, since the main goal of the present study is to in- 
vestigate the difference between TV-body and FP methods for ro- 
tating stellar systems, we need to apply the same criteria with the 
FP approach. The tidal boundary (or tidal energy) is adjusted dur- 
ing evolutio n as described in the previous studies (IKim et all2002l 
l2004ll2008l) . 

Table[2]shows the initial parameters used for the TV-body sim- 
ulations. For comparison, we also performed FP simulations of 
M2A models. There are several parameters that determine the clus- 
ter evolutio ns such as the concentration parameter Wo, the initial 
rotati on ujq ( Einsel & Spu rzemlll999l: IKim et alj|2002l; lErnst et al.l 
|2007|) . the mas s spectrum dKim et alj2004l;lKhalisi et alj200^ and 
tidal boundary dBaumgardd200lUKinj2003l ; lErnst et alj2007t) . We 
fixed Wo = 6, but varied ujo and mass spectrum to investigate the 
effect of the initial rotation and the mass spectrum on the cluster 
evolution. The initial rotation cjo are varied from 0.0 to 1.5. For the 
mass spectrum, we varied the individual mass ratio m^/mi from 
2 to 20 while the total mass ratio M\ /M 2 is fixed to 5 or 20. One 
of the most important processes of the stellar system with multiple 
mass components is the energy equipartition. The energy equipar- 
tition is a tendency for different mass components to have similar 
kinetic energies. However, in some cases, the equipartition is not 
completely ac hieved, which is called as the equipartition instability 
(Spitzeri [l969T) . To determine whether the equipartit ion happens to 
be achieved or not in two-component mass systems, ISpitzejd 1969b 
derived an analytic equipartition stability parameter, 



Mi 1 mi 



3/2 



(1) 



He suggested that the energy equipartition between low and high 
mass stars takes place when S < S C ru = 0.16. After Spitzer's 
study, many authors have studied the energy equipartition pro- 
cess of two-componen t systems by theoretical approaches (e.g., 
Light man & Fallll 1978h and by several numerical methods such as 
Mont e-Carlo approaches to solve the FP equation ( Spitzer & Hart 
1 197 ll) . the direct integratio n of FP equation dKim et al.l 1 199a) 
and TV-body simu lations dPortegies Zwa rt & McMillan 2000). 
IWatters etaUd2000l) performed Monte-Carlo simulations with var- 
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Table 1. Initial parameters for all models. 



Model 


W 


0)0 


r t id/r c 


T rot /\W\ 


7V*runs 


m 2 /mi 


Mi/M 2 


AT 2 


S 


A 






0.0 


18.0 


0.000 


20,000*3 
















0.3 


14.5 


0.035 


20,000*3 












M2At 
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0.6 


9.9 


0.101 


20,000*3 
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1818 


0.566 


1.056 
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n 1 
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1.2 


5.4 


0.196 


20,000*3 
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4.4 


0.222 


20,000*3 
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M2B 
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0.6 


9.9 


0.101 


20,000*1 
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769 


2.236 


9.518 






1.2 
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0.196 


20,000*1 
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20,000* 1 
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0.6 


9.9 


0.101 


20,000*1 


10 
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392 


6.325 


50.24 
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0.196 


20,000*1 
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18.0 


0.000 


20,000*1 












M2D 


6 


0.6 


9.9 


0.101 


20,000*1 


20 
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198 


17.89 


256.2 






1.2 


5.4 


0.196 


20,000*1 
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18.0 


0.000 


20,000*1 












M2Ae 


6 


0.6 


9.9 


0.101 


20,000*1 
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20 


488 


0.141 


0.264 






1.2 


5.4 


0.196 


20,000*1 













tFokker-Planck simulations are performed for comparison. 



ious two-component mass spectra and introduced an empirical 
equipartition stability parameter 



M2A 



A : 



Ma 
Mi 



m,2 
mi 



(2) 



and found that the critical value for energy equipartition is A cr u = 
0.32. Most of our models have S > S cr it and A > A cr i t . The ratio 
of the rotational kinetic energy to the potential energy is known 
to be a measure of the 'temperature' of the rotating system with 
higher value being called cold system. If this parameter is greater 
than 0.14, the system is known to become dyn amically unstable 
again st the formation of the bar-like structure (Ostri ker & Peebles! 
1 1973b . Thus, models with uiq greater than 0.9 are expected to evolve 
to elongated shape. We designate these models as rapidly rotating 
models and the other models with ojq ^ 0.6 as slowly rotating 
models in this study. 



3 SLOWLY AND RAPIDLY ROTATING CLUSTERS 

To investigate the effect of the initial rotation on the evolution of 
clusters, we focus on M2A models which have various amounts of 
the initial rotation. In Fig. [I] we show the time evolution of the ratio 
of the rotational kinetic energy to the potential energy, T rot /|W|. 
The dashed and solid lines represent FP and ./V-body results, re- 
spectively. The ini tial half-mass relaxati on time is measured as fol- 
lows suggested bv lSpitzer&Hartl ( ll97lh . 



r rh . = 0.138- 



(3) 



G 1 / 2 77T 1 / 2 lnA' 

where N,rh.o,G,m and InA are total number of stars, initial 
half-mass radius, gravitational constant, mean mass of stars and 
Coulomb logarithm, respectively. It is well known that a sys- 
tem with rigid-body rotation suffers a secular instability when 
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Figure 1. The ratio of rotational kinetic energy to potential energy. Solid 
and dashed lines represent A^-body and FP results, respectively. Upper 
panel shows the ratios for slowly rotating models and lower for those of 
the rapidly rotating models. The dotted lines denote a criterion of bar insta- 
bility lOstriker & Peeblesll 1973b . For rapidly rotating models the A-body 
and FP results show the significant difference due to the bar instability. 




T rot /\W\ is larger than 0.1 4 dOstriker & Peeblesl Il973l) . Later, 
from numerical simulations, ISellwoodl ( 119811) confirmed that the 
criterion is valid for more realistic rotation curves. Our models with 
initial value of T rot /\W\ < 0.14 are shown in the upper panel 
and those with T rot /\W\ > 0.14 are shown in the lower panel. 
The results of A^-body and FP show similar behaviors for the mod- 
els with initial Ti-ot/IVKj < 0.14. However, for the models with 
initial T rot /\W\ > 0.14, the initial evolution depends on differ- 
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Figure 2. Initial axis ratio of clusters with different initial rotations. Cross 
and plus symbols show intermediate (b/a) and minor (c/a) axis ratios, re- 
spectively. We estimate axis ratios with stars in half mass ellipsoidal radius. 
The clusters show more oblate shapes as the initial rotation increases. 



ent numerical approaches. The iV-body simulation represents more 
rapid decrease of T ro t/\W | with time than the FP simulation in the 
early phase. For the rapidly rotating models one can observe the 
construction of bar-like structure and the total rotational energy de- 
creases very quickly. Therefore, the FP approach seems to be not 
appropriate in describing the evolution of rotating models with ini- 
tial value of T ro t/\W\ > 0.14. 

In order to investigate the evolution of the cluster shape, we 
calculate axis ratios of clus ters by using the method suggested by 
iDubinski & Carlberddl99ll) . They defined a tensor, 

My = (4) 



with an ellipsoidal radius, 



<? = 



x 2 + 



+ 



(b/a)' (c/a) 



1/2 



(5) 



where a, b and c are axis lengths with a ^ b c. The axis ratios 
are derived from the tensor through 

1/2 / « ^ \ 1 / 2 

(6) 



M xx 



and — 



M z 



M x 



where M xx , M yy and M zz are the principal components of the ten- 
sor. In order to compute the tensor My , we need to know the axis 
ratios. Therefore, for the simultaneous determination of the tensor 
and axis ratios, we need to perform an iterative calculation. We, 
first, assume certain set of axis ratios (e.g., b/a = 1 and c/a = 1) 
and compute My, which gives another set of axis ratios. Obviously, 
the resulting axis ratios will be different from the assumed values, 
and therefore, can be used as an input for improved estimation of 
Mij . We carry out the iteration until the relative difference of axis 
ratios becomes less than certain criterion (we include a value of 
10 -4 for this study). Fig.|2]shows the axis ratios as a function of ini- 
tial rotation for M2A models. We calculate axis ratios with stars in 
the ellipsoidal radius including half-mass of the cluster. The shape 
of rotating cluster is oblate initially due to the initial rotation (i.e., 
b/a — 1 and c/a < 1). The minor axis ratio c/a decreases when 
the initial rotation increases from 0.3 to 1.5. 



Figure 3. Time evolution of intermediate b/a (black) and minor c/a (grey) 
axis ratios for M2A models. In the early stage for rapidly rotating models, 
the intermediate axis ratios decrease rapidly due to the bar instability. The 
shapes of clusters tend to become spherical after a few T r hfl due to the 
significant loss of angular momentum. 



In Fig. [3] we show the evolution of the axis ratios, b/a and 
c/a for M2A models with different initial rotations. The left pan- 
els are the result of slowly rotating models and right panels are for 
the rapidly rotating models. The minor axis ratios c/a increases 
with time for slowly rotating models because the cluster loses an- 
gular momentum. For the rapidly rotating models, the intermedi- 
ate axis ratio decreases at the beginning, because of the develop- 
ment of the bar instability. Due to this instability, cluster shapes 
become tri-axial or even prolate in a dynamical time scale which is 
much shorter than the relaxation time. The intermediate axis ratio 
also decreases during this phase. The decrease of both axis ratios is 
more rapid when the initial rotation becomes larger. When the ratio 
T rot /\W\ becomes smaller than 0.14 by the loss of the rotational 
energy as shown in Fig. [T] the bar instability disappears and the 
cluster becomes axisymmetric. 

Density contour maps on XY plane of clusters with different 
initial rotation parameters are shown in Fig. |4] The top and bottom 
panels show the initial shapes and the shapes at the core collapse, 
respectively. In the bottom panels, the size of contours show that 
the remaining mass of cluster at the core collapse becomes smaller 
as the initial rotation increases. Density structures for models with 
ujq = 0.9, 1.2 and 1.5 in middle panels represent shapes at the time 
of maximum elongation. One can clearly observe bar-like struc- 
tures of rapidly rotating models. Shapes become to be those of the 
axisymmertic systems at the time of core collapse as shown in Fig 
[3](i.e., b/a = 1 and c/a < 1). 



4 COMPARISON BETWEEN N-BODY AND FP RESULTS 
4.1 Mass and Energy 

In this section, we compare the evolution of ./V-body and FP simu- 
lations for M2A models with three different point of views: overall 
evolution (mass and energy), central evolution (central density and 
velocity dispersion) and rotational evolution. Fig.|5^ shows the evo- 
lutions of total mass. For the slowly rotating models, the time evo- 
lutions of total cluster mass for A^-body simulations agree well with 
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Figure 4. Projected density contour maps of M2A models with initial rotation u)q =0.0, 0.6, 0.9, 1.2 and 1.5 on XY plane. The top panels (a-e) show initial 
models and bottom panels (k-o) show density structures at the time of the core collapse. The middle panels of loq = 0.9, 1.2, 1.5 (h-j) are density structures 
at the time of maximum elongation. 
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Figure 5. Evolution of mass for M2A models: (a) Evolution of total mass. 
Dashed lines represent FP results, and solid lines show A-body results. The 
total mass decreases due to escaping stars through the tidal energy thresh- 
old. For rapidly rotating models, there are steep mass losses at the early 
stages because of the bar instability, (b) Time evolution of the mass compo- 
nents mi (red) and ni2 (blue). Total mass of m\ decreases faster than that 
of m,2- 



those of FP results. However, there exist significant differences be- 
tween FP and TV-body results for the rapidly rotating models. For 
an instance, TV-body calculations show significantly higher mass 
loss rate than FP results for rapidly rotating models with wo = 1.2 



and 1.5, especially in the very early times. The significant amount 
of mass loss at the very early stage also induces the large amount 
of angular momentum loss. After this stage, clusters have smaller 
number of stars with slower rotation and therefore evolve slowly 
compared to FP results. The time evolutions of total mass of each 
component are shown in Fig.|5p. Because the number of low mass 
stars is much larger than that of high mass stars, the evolution of 
total mass is similar to that of low mass stars. 



The total energy of the tidally limited stellar system decreases 
monotonically with time due to the escaping stars carrying ener- 
gies. Fig. [6^ shows the time evolution of the normalized total en- 
ergies. The results of TV-body and FP simulations agree well for 
slowly rotating models. While the total energy decreases slowly in 
pre-collapse stage, after core collapse the total energy decreases 
more rapidly due to the core expansion. For the rapidly rotating 
models, however FP and TV-body results show significant differ- 
ence because of the large amount of energy loss during the early 
stages for TV-body models. Fig. [6}} represents the time evolution 
of normalized specific energy (i.e., energy per unit mass) for each 
mass component. In overall, the specific energies of mi continue 
to decrease. On the other hand, the specific energy of m,2 increases 
slowly till the core collapse as a result of equipartition and mass 
segregation. However, the specific energies of 7712 decrease after 
the core collapse due to the mass loss through the tidal boundary 
because the process of mass segregation stops at the core collapse 
(i.e. the mean mass of a star in the central region increases until 
the core collapse and remains as a constant value after the core col- 
lapse; see §5.3). Note that the mean energies for both low mass and 
high mass stars immediately decrease for the model of uuq — 1.2 
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Figure 6. Upper panel shows the evolution of the total energy as a function 
of time. The diamond symbols represent the moment of the core collapse for 
Af-body results. Total energies decrease with time due to escaping stars. Af- 
ter the core collapse, energies decrease more rapidly becau se of core bounce 
due to the binary heating. In lower panel, there are the evolutions of the spe- 
cific energies for mi and ni2 . The specific energy of ni2 increases quickfy 
during the just prior of the core collapse because high mass stars move to 
the central region having deep potential due to the mass segregation. 
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Figure 7. Evolution of central density p c and velocity dispersion u c for 
M2A models. Dashed lines represent the FP results and the solid lines for 
the A^-body results. Different contrasts mean models with different initial 
rotations. p c and cr c of rapidly rotating models show significant differences 
between two approaches. 

and 1.5. Again, this is caused by the bar instability in the early 
phase of the evolution. 

4.2 Central Density, Velocity Dispersion and Core Collapse 

Fig. [7] shows time evolution of the central density, p c and the cen- 
tral velocity dispersion, a c obtained by using stars inside the core 
radius. The dashed and solid lines represent FP and TV-body results, 
respectively, and different contrasts represent models with different 



Table 2. Global evolution of M2A models for NBODY4 and mFOPAX re- 
sults 



Simul. 


W 


tcc/T r h.Q 


tdis/ T rh,a 




M cc 




0.0 


7.1 


22.0 


0.32 


0.80 




0.3 


6.7 


15.9 


0.42 


0.70 


mFOPAX 


0.6 


5.3 


9.2 


0.58 


0.52 




0.9 


3.8 


5.7 


0.67 


0.36 




1.2 


2.7 


3.7 


0.73 


0.23 




1.5 


2.1 


2.5 


0.84 


0.14 




0.0 


6.8 


21.2 


0.32 


0.80 




0.3 


6.2 


16.3 


0.38 


0.72 


NBODY4 


0.6 


5.7 


10.4 


0.55 


0.54 




0.9 


4.8 


7.1 


0.68 


0.35 




1.2 


3.6 


5.2 


0.69 


0.23 




1.5 


2.8 


3.7 


0.76 


0.15 



initial rotation. We confirm that the rot ation accelerates both the 
core collapse and cluster disruption (e.g JKim et al . 2008). In addi- 
tion, FP and TV-body results agree well for slowly rotating models. 
However, as the initial rotation increases the difference between 
FP and TV-body results becomes large due to the bar instability. 
Table [2] lists the core collapse times and the disruption times of 
M2A models. For example, the model with ujq = 1.5 reaches the 
core collapse at 2.\Trh,o and 3.7T r h,o for FP and A/-body calcula- 
tions, respectively. As the initial rotation uio increases, the cluster 
spends more time of its whole life in the pre-core collapse phase 
(i.e., t C c/tdia increases). For TV-body simulations, the ratio of the 
core collapse time to the disruption time t cc /tdi S of the model with 
u)q = 1.5 is ~3/4 while that of the model with wo = is ~l/3. 
We observe that the core collapse and disruption times of TV-body 
results are significantly longer than those of FP results for rapidly 
rotating models because these models are redefined as small and 
slowly rotating systems after the bar instability. Although there are 
significant differences between TV-body and FP time scales, the to- 
tal mass at the time of core collapse (designated as M cc ) of TV- 
body and FP results show a good agreement. The prese nce of the 
mass spectrum also accelerates the evolution of cluster dKim et all 
2004). We found that the core collapse times of M2A models are 
20—40% smaller than those of single mass systems considered by 
iKim et alj (2008J). On the other hand, the disruption times are only 
10% smaller. We also look into the evolution of clusters with var- 
ious mass spectra. Table [3] shows the collapse and disruption time 
scales for other TV-body models. As mijm\ increases, the evolu- 
tion of clusters is accelerated. From M2A to M2D models, the core 
collapse time decreases when 7772 /mi increases, but is less affected 
by the initial rotation. The disruption time, however, depends on 
both the mass spectrum and the initial rotation. 

Fig. [D shows the evolution of the central velocity dispersion 
(<r c ) of each mass component as a function of p c . a c of mi and 
7712 are divided into two parts and increase gradually until the core 
collapse. The total a c approaches that of 7772 because the fraction of 
high mass stars in the core increases with time and the core is finally 
filled with high mass starts due to the mass segregation. We can 
clearly see the equipartition as two distinct branches of a c versus 
p c . As a result of equipartition, a c of mi becomes about \/2-times 
greater than that of 7772 because we use individual mass ratio of 
t7i2 /mi = 2. Once the establishment of energy equipartition, the 
evolutions of a c of mi and 7712 are represented by simple power- 
law. During the post-core collapse stage, however, FP and TV-body 
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Table 3. Global evolution of other NBODY4 models 
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slow rotation M2A rapid rotation 



Model 


OJQ 




tdis/ T rh,0 




M cc 




0.0 


1.61 


10.18 


0.16 


0.94 


M2B 


0.6 


1.64 


5.01 


0.33 


0.79 




1.2 


1.38 


2.16 


0.64 


0.34 




0.0 


0.59 


5.92 


0.10 


0.97 


M2C 


0.6 


0.62 


2.79 


0.22 


0.89 




1.2 


0.70 


1.36 


0.51 


0.44 




0.0 


0.42 


3.59 


0.12 


0.96 


M2D 


0.6 


0.34 


1.65 


0.21 


0.91 




1.2 


0.36 


0.72 


0.5 


0.55 




0.0 


9.37 


27.56 


0.34 


0.73 


M2Ae 


0.6 


6.72 


12.86 


0.52 


0.49 




1.2 


4.05 


5.89 


0.69 


0.22 
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Figure 8. Central velocity dispersion cr c for mi (red), ni2 (blue) and all 
stars (black) as a function of central density p c . Dots show A-body results, 
and dashed lines mean FP results. Initially, cr c for both mass components 
are the same. They evolve into the different ways due to the equipartition. 
a c for low mass becomes about v^-times larger than that for high mass 
because of equipartition. 

results show some differences. iV-body results are more dispersive 
than FP results, unlike the early evolution. This is because there is 
the large fluctuation of p c as shown in Fig.[7J which is from the lack 
of stars within the core radius after core collapse. On the other hand, 
rapidly rotating models that lead to the bar instability have different 
evolutionary behaviors between FP and iV-body for the very early 
phase. Because the bar instability delays the relaxation processes, 
the mass segregation in TV-body takes place later compared to FP 
models. The gap between a c of mi and 7712 for iV-body results 
is slightly smaller than that of FP results. This means that the bar 
instability also affects the equipartition process (see §5.2 and Table 
|4]for more details). 

In Fig. [9] we show the time evolution of the Lagrangian radii 
for M2A models. We estimate Lagrangian radii containing 1,5, 10, 
20, 50 and 75 per cent of initial cluster mass. The dashed lines and 
the solid lines represent the FP and iV-body results, respectively. 
Although our models are flattened due to the initial rotation and 



Figure 9. Lagrangian radii including 1,5, 10, 20, 50, 75 per cent of initial 
cluster mass for M2A models. For rapidly rotating models, the results of TV- 
body and FP shows significant differences while those for slowly rotating 
models are similar. 

become triaxial by the bar instability, we use the spherical radius 
to estimate these Lagrangian radii for the simplicity of comparison 
between A^-body and FP simulations. We found that, for the slowly 
rotating models, the results of FP and A^-body show an excellent 
agreement except for the final stages when there are only small 
number of stars in the clusters. However, for the rapidly rotating 
models, the results of FP and A^-body are significantly different. 
Again, the bar instability and induced mass loss are the main reason 
for these differences. 



4.3 Rotational properties 

To understand the effects of initial rotations on the cluster dynam- 
ics, we investigate the evolution of angular momentum. Because we 
assume that stars escape through the tidal energy threshold, clusters 
lose their angular momentum continuously. Also, we expect ex- 
change of angular momentum between different mass components, 
similar to energy exchange. After encounters, stars that lose their 
energy spiral into the central parts while stars that gain the energy 
move outward. Also, loss of angular momentum makes the orbits 
of stars be eccentric but gaining angular momentum does the or- 
bits be less eccentric. Fig. [10] shows the time evolution of specific 
angular momentum which is defined as 

where fl, vi are relative positions and velocities of stars to the cen- 
ter of mass, respectively, and we integrate all stars in the cluster 
to compare A^-body and FP results. L z , spe c is z-direction compo- 
nent of L spec . Solid lines are the results of Af-body and dashed 
lines represent the results of FP. Red and blue colors show low and 
high mass stars, respectively. For slowly rotating models, the re- 
sults of A^-body and FP are similar. The specific angular momenta 
decrease monotonically. Due to the mass segregation, high mass 
stars that encounter with many low mass stars migrate to the inner 
part where rotational velocity is smaller than that of the outer part 
of the cluster (see rotation curves in Fig. [I2ji and Fig. HBll. This 
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Figure 10. Time evolution of the specific angular momentum L ZiSpec . 
Solid lines are the results of FP and dashed lines are the results of N-body. 
Red and blue colors mean low mass and high mass stars, respectively. Di- 
amond symbols represent the moment of the core collapse. For slowly ro- 
tating models (left panels), L z , apec of mi and rrt2 are divided. However, 
Lz,spec of mi and m,2 decrease quickly compared to the results of FP for 
rapidly rotating models (right panels). 
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Figure 11. Distribution of the rotational velocities of stars as a function of 
cylindrical radius. The cluster's size and degree of the anisotropy decrease 
with time. 



is also related to the angular momentum exchange between differ- 
ent mass components (see §5.4 for more details). Thus, the specific 
angular momentum of m,2 becomes smaller than that of mi during 
the evolution. The difference of angular momentum between two 
mass components is most prominent at the core collapse, where the 
core collapse time is marked as diamond symbols in Fig. [10] For 
the models with ui — 1.2 and 1.5, the angular momentum de- 
creases rapidly during the early stage. In addition, there is no clear 
split of the specific angular momentum of low and high mass stars 
due to the rapid evolution induced by the bar instability. Interest- 
ingly, the model with ljq = 0.9 shows the transitional evolution. 
While the model still has a triaxial shape (T < 1.5r r h,o, see Fig. 
[3J, the evolution is similar to that of the models with uio = 1.2 and 
1.5. However, after the shape becomes the axisymmetric again, the 
evolution is similar to that of the models with Uo = 0.3 and 0.6. 

4.3.1 Slow rotation (uo = 0.6) 

To investigate the evolution of rotational properties of our cluster 
models, we focus on rotational velocities, V$ of two models in de- 
tail. As a representative of slowly rotating models, we chose a ro- 
tating model with ujq = 0.6 of M2A. Fig.[TT]shows the distribution 
of V(j> at four different epochs T = 0,3, 5.7 and 8r r h,o- The core 
collapse time is T = 5.7T r h.o- We combined results of three runs 
with N =20,000. At T = 0, the distribution is asymmetric to the 
positive direction due to the initial rotation. The dispersion of 
is large at the center and becomes smaller along the cylindrical ra- 
dius, R = (x 2 + y 2 ) 1 / 2 . As the cluster loses its mass and angular 
momentum, the size of cluster and the degree of asymmetry in V4, 
decrease with time. Finally, the distribution of becomes sym- 
metric compared to the initial distribution. The rotation curves of 
both mass components are shown in Fig. I12fa . The rotation curves 
of mi and m,2 are identical at T = because we assume that 
the distribution of positions and velocities for both mass compo- 
nents are the same initially. At T = 37VA,o> the rotation velocity of 
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Figure 12. Upper panel shows rotation curves at time T = 0,3,5.7 
(core collapse) and 8T r ^ q. Solid, dotted and dashed lines mean the rota- 
tion curves of all stars, high mass stars and low mass stars, respectively. At 
time T = 0, rotation curves are similar between low and high mass com- 
ponents. However, after few t,/i 0, the rotation curve of mi drops while 
that of rrt2 remains. The curve of z-direction angular momentum is shown 
in lower panel. Angular momenta also decrease with time. 



mi becomes smaller at all radii because low mass stars with large 
angular momenta escape from the system. On the other hand, the 
rotation curve of m,2 remains as that of the initial curve for longer 
time. At T = T cc , the rotation velocity of m,2 also decreases but 
is still higher than that of mi. The cluster rotates very slowly at 
T = 8r r h,o- At this time, the rotation curves of mi and m,2 are 
flattened compared to other epochs. It is interesting to note that the 
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Figure 13. Time evolution of z-direction specific angular momentum 
Lz,spec (upper), L Z: spec of inner half-mass (middle) and L ZtSpec of outer 
half-mass (lower). Though the angular momentum curve of high mass com- 
ponents is larger than low mass components at all radii, L ZiSpec of high 
mass components in whole cluster becomes smaller than low mass compo- 
nents after few r r j, q. This is because most of high mass components are 
located in central region which rotates slowly due to the mass segregation. 
The evolution of the inner region shows similar behavior with that of whole 
cluster. On the other hand, for the outer region, L z , spec of high mass com- 
ponents is higher than that of low mass components always. 



peak position remains at a constant value measured in the units of 
half-mass radius while the peak rotation velocity decreases. 

Fig. 1 12b shows the radial profiles of mean angular momentum 
for mi and m z . The angular momentum curves of mi and 7712 are 
identical at T = and L z of mi decreases more rapidly com- 
pared to that of 77i2. Initially, the curve is a power-law with index 
of 1.5 within (e.g., the power law index is 2 for a rigid body ro- 
tation.). The power law index becomes smaller and goes to one at 
T = ST r h,o because the rotation curve becomes flatter. In Fig. 1131 
we show evolutions of specific angular momentum L Z:Spec in dif- 
ferent radial ranges: in the whole cluster, within half-mass radius 
and outside of half-mass radius. L ZtSpec decreases with time due 
to the escaping stars with angular momenta. In the entire cluster, 
L Z:Sp ec of mi and 7712 decrease together until T — 2r r h,o- Al- 
though L z of 771,2 is larger than that of mi at the all radii as shown 
in Fig. 112b. L Z:Spec of 7712 becomes smaller than that of mi be- 
cause high mass component tends to be more concentrated in the 
central region than low mass stars due to mass segregation (see §5.3 
for more details). Within the half-mass radius, L ZtSpec of mi de- 
creases continuously while that of 777,2 remains at a constant value 
until T ~ 0.5r r h,o because the rotation curve of high mass stars 
remains as the initial curve for a while as shown in Fig. I12fc and 
Fig. 112b. After a few r r h,o, the evolution of inner region follows 
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Figure 14. Distribution of the rotational velocities of stars as a function of 
cylindrical radius for the model with ujq = 1.5. Initially, most of stars have 
positive rotational velocities. However, the asymmetry rapidly decreases 
within airtj. 

that of entire cluster due to the mass segregation. For the outer re- 
gion, L Z:Spec of 77i2 is higher than that of mi throughout the whole 
evolutionary phase. Even for the mass segregation, the mean mass 
of outer region (see Fig|19) does not change significantly. It means 
that the mass fraction of the outer region is less affected than that of 
inner region. Thus, the mean rotation of 771,2 is still faster than the 
rotation of mi beyond half-mass radius. Finally, the cluster loses 
most of its mass and angular momentum at the end of evolution. 

4.3.2 Fast rotation (wo = 1.5) 

We also investigate the evolution of rotational properties for the 
model with uio = 1.5. As we mentioned earlier, this model is un- 
stable against bar instability and the shape of cluster quickly be- 
comes a prolate with the largest elongation at T ~ O.Itv^o- In Fig. 
[141 we show the distribution of tangential velocities of stars at 
T = 0,0.1, 1 and 2.8r r h,o (core collapse). Initially, the distribu- 
tion is more skewed toward the positive direction than the model 
with ojo — 0.6 as shown in Fig.QT] Only less than 10% of stars 
have negative value at T — 0. The distribution becomes symmet- 
ric and also the size of cluster becomes smaller as similar to the 
model with ujq — 0.6. At core collapse, only about 10% of stars re- 
main in the cluster and the cluster rotate slowly. Fig. [Tsb shows the 
rotation curve at T — 0, 0.1, 1 and 2.8r r h,o- The rotation curves 
of mi and 7712 decrease together, indicating that the bar instabil- 
ity (i.e., large mass, energy and angular momentum loss) disturbs 
the relaxation processes. The peak of rotational velocity decreases 
rapidly with time. However, the peak position at T = O.Itv^o is 
slightly larger than the other epochs due to the effect of bar from 
the instability. At core collapse, the rotation curve is nearly flat but 
still remains, though the cluster lose most of the mass and the an- 
gular momentum. The mean angular momentum along the radius 
is shown in Fig. 1 15b. Similar to the result of model with ojo = 0.6, 
it shows a power-law distribution. Initially, the power law index is 
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Figure 15. Rotation curves at time T = 0, 0.1, 1 and 2.8 (core collapse) 
Trh,o with u>o = 1.5 (upper). Lines have same meaning to those of Fig. 
1121 At time T = 0, rotation curves are similar between low and high mass 
components. Unlike rotation curves of in Fig |12l those of low and high mass 
components are not divided much. Curves of z-direction angular momen- 
tum are shown in lower panel. 
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Figure 17. Maximum mass evaporation rate £ e ,post after core collapse as 
a function of the individual mass ratio. Error bars show standard deviations 
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Figure 16. Dimensionless mass evaporation rate. X-axis means the evolu- 
tionary phase (i.e.; pre- or post-core collapse). For uo = 0.0 and 0.6, the 
result of N -body and FP are well agreed. On the other hand, for ojq = 1.2, 
there is a spike meaning large mass loss induced by the bar instability in 
A-body result. 



~2 within like a rigid body rotation but it becomes close to 1 at 
core collapse. 



5 DISCUSSION 

5.1 Mass evaporation 

To investigate the evolution of mass in detail, we first define the 
dimensionless mass evaporation rate such as 

T rh dM 
M dt 



of data of , 



\e,post 



for different mass spectra. 



(8) 



where r r h and M are the half mass relaxation time and total mass 
of a cluster at time T, respectively. Fig. [T6] shows the behavior of 
£ e of TV-body and FP simulations for M2A models with ujq =0.0, 
0.6 and 1.2. We divide the evolution into pre- and post-core col- 
lapse phases to investigate the evolution of £ e more clearly. The 
mass evap oration rates are known to be cons tant for self-similar 
case (e.g., lHenorJll96ll : llxe & Ostrikerlll987t) . However, the rate 
changes with time because our models are not self-similar. In early 
phase, £ e increases with the initial rotation. For slowly rotating 
models, we see a very good agreement between A-body and FP 
results in pre-core collapse phase. On the other hand, for model 
with ujq = 1.2, there is a significant difference between results of 
A^-body and FP. The A^-body results show a spike at the early time 
while FP results decrease monotonically. This spike is related to 
the rapid mass loss induced by the bar instability. After the spike, 
£ e of A^-body suddenly decreases below £ e of FP result. It is in- 
teresting that £ e ,cc, the mass evaporation rate at the core collapse, 
of A^-body and FP show similar results even though there is a big 
difference before core collapse. After core collapse, £ e of A^-body 
and FP increase toward a peak value and decrease afterward. The 
large differences of £ e between TV-body and FP results at the end 
of the evolution is due to small num ber of rema i ning s tars and thus 
do not have statistical significance. iKim et all 120021 ) carried out 
FP simulations for rotating clusters with single mass system. They 
also calculated £ e with different initial rotation ujo =0.0, 0.3 and 
0.6. We notice that the evolutionary shapes are similar between 
single mass and the 2-component mass systems. However, £ e for 2- 
component mass systems show about 30% enhancement compared 
to £ e for single mass systems in pre-core collapse phase. This en- 
hancement could have been induced by the ener gy exchange pro- 
cess in multi-component models, as noticed by iLee & Goodman! 

They calculated f e with various initial mass functions and 
found that £ e increases when the cut-off mass ratio (i.e., m,f/mi 
if the mass ranges from m, to mf) increases. To confirm the re- 
lationship between the mass evaporation rate and the mass ratio, 
we compute the maximum evaporation rate after the core collapse 
£e,post for A^-body results with various mass spectra. Fig.ll7lshows 
£e, P ost as a function of the individual mass ratio m,2 /mi . As shown 
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Figure 18. Examples for the evolution of the equipartition parameter § e 



Dashed lines mean the minimum values £ eq 
core collapse time. The determination of £ e 



, , and dotted lines show the 
.„ is mentioned in text. 



in Fig. [T6j the peak £ e in post-core collapse phase are very sim- 
ilar with different initial rotations, so we average the results of 
each mass spectrum for different initial rotations. £ e ,p OS t of M2A- 
D models increases with increasing 7712/mi and follows a simple 
power-law. For M2Ae model, £, e , P ost is slightly smaller than that 
of M2A model because the fraction of high mass stars for M2Ae 
model is small. 



5.2 Energy equipartition 

As presented in §4.2, high and low mass stars in the core approach 
to the 'thermal' equilibrium state by the two body relaxation. To 
investigate the energy equipartition in detail, we adopt the equipar- 
tition parameter 



m 2 a 2 
raicr, 2 



(9) 



like previous studies dWatters et al 1 120001: iKhalisi et alJl2007t) . We 
calculated £ eq for stars inside the core radius. Fig. II 8 1 shows the 
evolution of £ e<3 for models M2Ae and M2D without initial rotation 
for examples. Dotted lines represent the core collapse time. For 
the M2Ae model, £ eq approaches to unity and becomes less than 
unity at some moments. However, for the M2D model, £ e9 never 
approaches to unity. The data is very noisy after the core collapse 
because only a few low mass stars remain in the core. To determine 
the minimum value £eq,min, we find the best polynomial fitting 
function by varying the order from 5 to 15. 

Table [4] lists (, eq . m in of all models. Only M2Ae models with 
wo = 0.0 and 1.2 have lower ( e ,, m i„ than 1.05 which is the value 
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Figure 19. Mean mass in mass shells between different Lagrangian radii 
with different initial rotations. Mean masses of inner shells increase with 
time before the core collapse while mean masses of outer shells decrease. 
This shows the process of the mass segregation. After the core collapse, 
mean masses of other shells also increase with somewhat time gaps but this 
is not from the mass segregation but because low mass stars escape more 
frequently than high mass stars. 



used in IKhalisi et alj 120071) for the energy equipartition. On the 
other hand, £ eq ,min increases when the equipartition instability pa- 
rameters S or A become larger. f et j,mi« for slowly rotating mod- 
els (i.e., ljo ^ 0.6) are similar to each other. There are, however, 
significant differences of £ C9 ,min between slowly rotating models 
and rapidly rotating models. For M2A models with uio 0.9, 
€,eq,min increases gradually with initial rotation. This is another 
phenomenon of the bar instability obstructing the relaxation pro- 
cess. 

5.3 Mass segregation 

A simple met hod measuring th e degre e of the mass segregation is 
suggested by|G iersz & Heggie ( 1996). They calculated the mean 
mass in a space between different Lagrangian radii (i.e., Lagrangian 
shells) to see the chan ge of mass d i strib ution in un-equal mass sys- 
tems. More recently. IKhalisi et al.1 1200% carried out iV-body sim- 
ulations with different mass spectra and found that the mass segre- 
gation occurs inward direction (i.e., the mean mass of each shell is 
decoupled stepwise from outside, see models A and B in Fig. 6 of 
IKhalisi et all (20071) ). Fig.[l9]shows the evolution of mean mass in 
different Lagrangian shells as a function of time. In pre-core col- 
lapse phase, due to the mass segregation, the mean mass of the in- 
nermost shell increases while mean masses decrease in outer shells. 
Note that the maximum mean mass of the innermost shell does not 
depend on the initial rotation and also the innermost shell is nearly 
fully-occupied by high mass stars (i.e., < m >~ m 2 — 1.833) 
at the time of core collapse. The mean mass of the innermost 
s hell stays at a constant v alue after the core collapse. According 
to[G iersz &Heggid lll996), after the core collapse, mean masses of 
inner shells with r < r 2 o% slightly decrease because high mass 
stars are removed by binary formation and the mass distribution fi- 
nally reaches a steady state in Lagrangian coordinate although the 
system expands with time. Our simulations, however, show s ignifi- 
cant differences from the previous studies (IGiersz & Heggie 1996; 
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Table 4. f e 



, of all models 
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S 


A 






Model 


S 


A 


wo 




M2A 


0.566 


1.056 


0.0 


1.105 


M2C 


6.325 


50.24 


0.0 


3.017 
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1.091 








0.6 


2.923 








0.6 


1.107 








1.2 


3.314 
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1.158 








0.6 


5.224 








1.5 


1.201 








1.2 


6.672 
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Figure 20. Evolution of the angular momentum for the model with uo = 
0.6. When we consider the angular momentum of escapers, the total angular 
momentum is conserved (thick solid line). The total angular momentum of 
mi slightly increases (thick dotted line), while that of ni2 decreases (thick 
dashed line). This shows the existence of the angular momentum exchange 
from m,2 to mi. The words 'esc' and 'exc' mean the degree of the angular 
momentum loss by escape and exchange, respectively. 



Figure 21. Detailed evolution of the angular momentum of rri2 including 
those of escapers (i.e., thick dashed line in Fig. 1201 as a function of t norm . 
tnorm is a time unit normalized by half-mass relaxation time (see the text). 
Open circles represent the iV-body results. At the end of evolution, one can 
distinguish the amount of angular momentum loss by exchange and escape. 
The solid line shows the Eq. (11) with £ e xc in TablefS] Note that the line is 
shifted as much as the amount of angular momentum loss by escape. 



iKhalisi et al]|2007t) in post-core collapse phase. In Fig. [19] mean 
masses of inner shells continue to increase even after the core col- 
lapse. This is because our models are tidally-limited. Low mass 
stars escape from the cluster more rapidly than high mass stars as 
shown in Fig. [5] As most of low mass stars escape, mean masses of 
outer shells increase at the end of evolution. 

5.4 Angular Momentum Exchange 

Because our TV-body model includes two different mass compo- 
nents, we are able to investigate the angular momentum exchange 
between different mass components. This has not been studied 
carefully yet and is an important subject to understand the evolu- 
tion of rotating star clusters. However, it is not easy to analyze the 
angular momentum transfer with the tidal boundary because total 
angular momentum of cluster decreases continuously by escaping 
stars. To distinguish the loss of angular momentum of a cluster be- 
tween escaping and exchange, we register the positions and veloc- 
ities of each escaping stars at the moment of escape. Fig.l20lshows 
the time evolution of the total angular momentum for a model with 



ujo = 0.6. Thin lines mean the total angular momentum of stars 
within the cluster and thick lines mean those of stars including es- 
capers. From the thick solid line, we observe that the total angular 
momentum including escapers is conserved as expected. Interest- 
ingly, the angular momentum of m,2 including escapers decreases 
while that of mi increases. Therefore, we conclude that there is a 
transfer of angular momentum from rri2 to mi. 

If the exchange of angular momentum is due to the two-body 
relaxation, we can define the angular momentum exchange rate as 
follows 

C _ _ T rh (t) dLz,m 2 -fm 1 , . 

Seic — ~f "77 ) \ l,J ) 

where L z , m2 _> mi is the amount of remaining angular momentum 
expected to be transferred from to m\. If the angular momen- 
tum transfer rate is a constant, the above equation can be integrated 
to give, 

(0)e"«— Wm (11) 

when we use a new time unit normalized by half-mass relaxation 
time tnorm = J 1/ T r h(t)dt. This rate can be a good parameter 
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Table 5. Angular momentum loss by escape and exchange for M2A models 



^0 L Zi o(m2) escape (%) exchange (%) £ exc 



0.3 


0.0283 


57.8 


42.2 


0.17 


0.6 


0.0448 


66.2 


33.8 


0.19 


0.9 


0.0555 


70.4 


29.6 


0.23 


1.2 


0.0571 


88.1 


11.9 


0.34 


1.5 


0.0635 


91.3 


8.7 


0.71 



to measure the efficiency of the angular momentum transfer be- 
tween different mass components. The angular momentum of clus- 
ters goes to when clusters are disrupted. So we divide the loss of 
angular momentum of mi into two contributions: through escap- 
ing and through transferring to mi. We also estimate the fractional 
angular momentum loss by two different processes. In Fig. |2T] the 
detailed evolution of the angular momentum of including those 
of escapers is represented by open circles. At the end of evolution, 
whole amount of angular momentum of 7712 disappears by escape 
or exchange. The solid line which is from Eq. (11) with suitable 
value of £ exc and shifted as much as the amount of angular mo- 
mentum loss by escape agrees well with the A^-body result. Table 
[5] shows the initial total angular momentum of 7712, the fraction of 
angular momentum loss by escaping and exchanging and the an- 
gular momentum exchange rate. Although the angular momentum 
exchange rate increases with the initial rotation, the fraction of an- 
gular momentum loss by exchange decreases. For the model with 
ljo — 1.5, the fraction is less than 10%. Rapidly rotating models 
have larger angular momentum exchange rate than slowly rotating 
models, but their lifetimes are very short compared to slowly rotat- 
ing models. Therefore, these rapidly rotating models do not have 
enough time to exchange angular momentum from 7712 to mi. 



6 SUMMARY AND CONCLUSION 

We performed numerical simulations of rotating stellar system with 
two mass components using NBODY4 and mFOPAX codes. By 
considering various mass spectra, we confirmed that both the initial 
rotation and the mass spectrum accelerate the evolution of the stel- 
lar system, as presented prev ious studies dEinsel & Spurzemll 19991 ; 
iKim et alj|2002 . |2004 I2008T) . However, we found that the initial 
rotation does not affect the evolution before the core collapse when 
the individual mass ratio 777,2 /mi is large enough. The mass evapo- 
ration rate is closely related to the acceleration of the evolution and 
increases with 7712 /mi . 

According to the instability criteria from lOstriker & Peeblei 

dl973h . we classified our models to slowly rotating mod- 
els (i.e., T ro t/\W\ <0.14) and rapidly rotating models (i.e., 
Trot /|W^| >0.14). By comparing the results of different ap- 
proaches, A^-body and FP simulations, we confirmed that two ap- 
proaches agree well with small differences on the time scales for 
the slowly rotating models. On the other hand, for rapidly rotating 
models, there are significant discrepancies between A^-body and 
FP results. From the investigation of shape of systems, we revealed 
that the bar instability happens at the beginning for rapidly rotating 
models in A^-body simulations. This bar instability induces unex- 
pected phenomena like the rapid loss of mass, energy and angular 
momentum. In addition, the bar instability hinders the two-body 
relaxation process, so the dynamical evolution of rapidly rotating 
systems is delayed as compared with FP results. We therefore con- 



cluded that the 2 dimensional FP approach is not valid for rapidly 
rotating cases because 2 dimensional FP approaches are unable to 
treat non-axisymmetric models. 

As the result of two-body interactions, low and high mass 
stars exchange their kinetic energies and happen to have similar 
kinetic energies. We confirmed t hat our models agree well with the 
equipartition instability criteria dSpitzerlll969l;IWatters et al.ll2000h 
for slowly rotating models. When the mass ratio becomes larger, 
it is hard to reach the complete equipartition state. Moreover, the 
equipartiton process is more disturbed for rapidly rotating models 
which suffer the bar instability. We also observed the exchange of 
angular momentum between low and high mass stars by investi- 
gating escapers and defined the angular momentum exchange rate 
Cezc- ^eic increases when the initial rotation increases. However, 
the amount of transferred angular momentum from high mass stars 
to low mass stars decrease because clusters with rapid initial ro- 
tation survive rather shortly compared to those with slow initial 
rotation. 
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